Homework 2

Exercise 1: ggplot

The following url: “http://cdiac.ess-dive.lbl.gov/ftp/ndp030/CSV-FILES/nation.1751_2014.csv” contains data on fossil fuel emissions.

1. Read data to R. Note that rows 1-3 contain information on the dataset itself. Delete these rows as they do not contain relevant information.

emissions <- read.csv("http://cdiac.ess-dive.lbl.gov/ftp/ndp030/CSV-FILES/nation.1751_2014.csv")
emissions <- emissions[4:nrow(emissions), ]
head(emissions)

2. Compute the total yearly \(CO_2\) emissions (column “Total.CO2.emissions.from.fossil.fuels.and.cement.production..thousand.metric.tons.of.C.”) summed over all countries (the world total \(CO_2\) emission). You can use a forloop over “years” or dplyr functions.

library(dplyr)
emissionsTotal <- emissions %>% 
  group_by(Year) %>%
  summarise(
    world_total = 
      sum(Total.CO2.emissions.from.fossil.fuels.and.cement.production..thousand.metric.tons.of.C.))

3. Plot the world (summed over all countries) \(CO_2\) emission over time in billion tonnes (Gt) per year, i.e. divide the quantity computed in (b) by 10^6. You can use a line or a scatter plot.

library(ggplot2)
ggplot(emissionsTotal, aes(x = Year, y = world_total/10^6)) +
  geom_line() + theme_bw()

4. Now read the dataset located at “https://raw.githubusercontent.com/lukes/ISO-3166-Countries-with-Regional-Codes/master/all/all.csv” which contains an assignment of countries to regions. Merge emissions dataset to the countries dataset. Note that in emissions dataset countries are given with all caps, but not in countries dataset. You need to change that before merging the two tables. Hint: use the function toupper(). Add a column ‘co2_emission’ equal to \(CO_2\) emission in Gt, i.e. ‘Total.CO2.emissions.from.fossil.fuels.and.cement.production..thousand.metric.tons.of.C.’/10^6

countries <- read.csv("https://raw.githubusercontent.com/lukes/ISO-3166-Countries-with-Regional-Codes/master/all/all.csv")
countries$name <- toupper(countries$name)
head(countries)

df <- emissions %>% 
  left_join(countries, by = c("Nation" = "name")) %>%
  mutate(co2_emission = Total.CO2.emissions.from.fossil.fuels.and.cement.production..thousand.metric.tons.of.C./10^6)

Column `Nation`/`name` joining factor and character vector, coercing into character vector

head(df)

5. Use dplyr to compute total annual \(CO_2\) (‘co2_emission’) emission per ‘sub.region’.

df <- df %>%
  group_by(Year, sub.region) %>%
  summarise(co2_emission = sum(co2_emission))
head(df)

6. Use ggplot to generate a stacked density plot the annual \(CO_2\) (in giga tonnes) by world regions (“sub.region”). Your plot should resemble something like this (but with other regional categories, and slightly different values). Hint: use geom_area() function with suitable parameters.

Source: https://ourworldindata.org/grapher/annual-co-emissions-by-region

ggplot(df, aes(fill = sub.region, x = Year, y = co2_emission)) +
  geom_area(stat = "identity") +
  theme(legend.position = "bottom")

Exercise 2: Gene expression data

In this exercise we will use the DNA microarray gene expression data. You can read more about it on page 5 of “The Elements of Statistical Learning”.

microarray <- read.table("https://web.stanford.edu/~hastie/ElemStatLearn/datasets/nci.data")
info <- read.table("https://web.stanford.edu/~hastie/ElemStatLearn/datasets/nci.info.txt",
                   skip = 12)
colnames(microarray) <- info$V1
microarray <- as.matrix(microarray)
microarray[1:5, 1:5]

        CNS      CNS    CNS  RENAL BREAST
[1,]  0.300 0.679961  0.940  0.280  0.485
[2,]  1.180 1.289961 -0.040 -0.310 -0.465
[3,]  0.550 0.169961 -0.170  0.680  0.395
[4,]  1.140 0.379961 -0.040 -0.810  0.905
[5,] -0.265 0.464961 -0.605  0.625  0.200

In the ‘microarray’ matrix columns correspond to samples, and rows to genes.

1. Randomly subset microarray to 500 genes. Then, plot a heatmap without clustering/dendrograms. You can use ‘asp = 0.2’ argument to change the aspect ratio of the heatmap.

set.seed(123456)
idx <- sample(1:500)
heatmap(microarray[idx, ], Rowv = NA, Colv = NA, asp=0.2)

2. Plot the previous heatmap with red/green color scheme. For your convenience here is the color vector you might like to use:

redgreen <- c("#FF0000", "#DB0000", "#B60000", "#920000", "#6D0000",
              "#490000", "#240000", "#000000", "#002400", "#004900",
              "#006D00", "#009200", "#00B600", "#00DB00", "#00FF00")

heatmap(microarray[idx, ], Rowv = NA, Colv = NA, asp=0.2, col = redgreen)

Now, plot the same graph but with dendrogram for rows (rows clustering).

heatmap(microarray[idx, ], Colv = NA, asp=0.2, col = redgreen)

4. Now instead of base heatmap, use interactive heatmaply() to generate the previous plot. You might want to add a command similar to the following %>% layout(margin = list(l = 150, b = 350), autosize = F, width = 600, height = 800) to the plot to set margins and to resize it.

library(heatmaply)
heatmaply(data.frame(microarray[idx, ]), Colv = FALSE, colors = redgreen) %>%
  layout(margin = list(l = 150, b = 350), autosize = F, width = 600, height = 800)

5. What interesting patterns do you observe? Are there some differences between conditions? Are some genes up down regulated for certain groups? No need for long answers just look at the heatmap state what you see.

Exercise 3: Hypothesis testing

Recall the movies data-frame we used in for lecture 3 exercises. It contains information on movies from the last three decates, which was scrapped from the IMDB database.

url <- "https://raw.githubusercontent.com/Juanets/movie-stats/master/movies.csv"
movies <- tbl_df(read.csv(url))
movies

1. Generate a boxplot of runtimes for action movies and commedies with jittered points overlaid on top. You might consider setting collor, fill and alpha arguments to modify clarity and transparency of the plot.

ggplot(movies %>% filter(genre %in% c("Action", "Comedy")), 
       aes(x = genre, y = runtime)) +
  geom_boxplot() +
  geom_jitter(height = 0,alpha = 0.2, color = "grey30") 

2. Test a hypothesis that the action movies have higher mean runtime (length) than the comedies. Is the difference statistically greater than zero at significance level \(\alpha = 0.05\)?

t.test(formula = runtime ~ genre,
       data = movies %>% 
         filter(genre %in% c("Action", "Comedy")), 
       alternative = "greater")

Exercise 4: linear model

1. Read the data from “http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv” containing information on sales of a product and the amount spent on advertising using different media channels.

2. Generate a scatterplot of sales against the amount of TV advertising and add a linear fit line.

3. Now make a 3D scatterplot with axes corresponding to ‘sales’, ‘TV’ and ‘radio’.

4. The dataset has 200 rows. Divide it into a train set with 150 observations and a test set with 50 observations, i.e. use sample(1:200, n = 150) to randomly choose row indices of the advertising dataset to include in the train set. The remaining indices should be used for the test set. Remember to choose and set the seed for randomization!

set.seed(12345)
idx <- sample(1:nrow(advertising), 150)
train <- advertising[idx, ]
test <- advertising[-idx, ]

5. Fit a linear model to the training set, where the sales values are predicted by the amount of TV advertising. Print the summary of the fitted model. Then, predict the sales values for the test set and evaluate the test model accuracy in terms of mean squared error (MSE).

fit1 <- lm(sales ~ TV, data  = train)
summary(fit1)

Call:
lm(formula = sales ~ TV, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.5442 -1.9976 -0.0122  1.9375  7.4951 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7.236487   0.529751   13.66   <2e-16 ***
TV          0.045090   0.003164   14.25   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.266 on 148 degrees of freedom
Multiple R-squared:  0.5785,    Adjusted R-squared:  0.5756 
F-statistic: 203.1 on 1 and 148 DF,  p-value: < 2.2e-16

yhat <- predict(fit1, test)
(mse1 <- mean((test$sales - yhat)^2))

[1] 10.75636

6. Fit a multiple linerar regression model including all the variables ‘TV’, ‘radio’, ‘newspaper’ to model the ‘sales’ in the training set. Then, compute the predicted sales for the test set with the new model and evalued the MSE.
   Did the error decrease from the one correspodning to the previous model?

fit2 <- lm(sales ~ TV + radio + newspaper, data  = train)
summary(fit2)

Call:
lm(formula = sales ~ TV + radio + newspaper, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.7941 -0.9224  0.2934  1.1735  2.5006 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.796684   0.365948   7.642 2.61e-12 ***
TV          0.045314   0.001624  27.908  < 2e-16 ***
radio       0.190504   0.009962  19.124  < 2e-16 ***
newspaper   0.002500   0.006671   0.375    0.708    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.676 on 146 degrees of freedom
Multiple R-squared:  0.8905,    Adjusted R-squared:  0.8883 
F-statistic: 395.9 on 3 and 146 DF,  p-value: < 2.2e-16

yhat <- predict(fit2, test)
(mse2 <- mean((test$sales - yhat)^2))

[1] 2.98811

7. Look at the summary output for the multiple regression model and note which of the coefficient in the model is significant. Are all of them significant? If not refit the model including only the features found significant. Which of the models should you choose?

fit3 <- lm(sales ~ TV + radio, data  = train)
summary(fit3)

Call:
lm(formula = sales ~ TV + radio, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8650 -0.9098  0.2910  1.1932  2.5064 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.840365   0.345867   8.212 1.02e-13 ***
TV          0.045326   0.001619  28.002  < 2e-16 ***
radio       0.191739   0.009373  20.457  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.671 on 147 degrees of freedom
Multiple R-squared:  0.8904,    Adjusted R-squared:  0.8889 
F-statistic: 597.3 on 2 and 147 DF,  p-value: < 2.2e-16

yhat <- predict(fit3, test)
(mse3 <- mean((test$sales - yhat)^2))

[1] 2.962326

Exercise 4: classification

Email spam ESL

Exercise 5: clustering

Exercise 6: ggmap

1. Plot housing market prices
2. Plot route to work

Fire smoke https://cran.r-project.org/web/packages/PWFSLSmoke/vignettes/Maps_and_Timeseries_Plots.html Fire smoke https://cran.r-project.org/web/packages/PWFSLSmoke/vignettes/Data_Model.html

https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/ggmap/ggmapCheatsheet.pdf

http://web.stanford.edu/~imalone/VAM/ggmap.html

http://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html

http://urbanspatialanalysis.com/dataviz-tutorial-mapping-san-francisco-home-prices-using-r/

https://takeout.google.com/settings/takeout

https://medium.com/@tejasrr19/visualize-your-google-location-history-a2343b14a6fe https://shiring.github.io/maps/2016/12/30/Standortverlauf_post

https://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/heatmap/